Tarea 4
Tarea 4: Bayesian Inference Part I
Integrantes :
- Alumno 1: Camila Bergasa
- Alumno 2: Matías Seda
Cuerpo Docente:
- Profesor: Felipe Bravo M.
- Auxiliar: Sebastian Bustos e Ignacio Meza D.
Fecha límite de entrega:
Índice:
Objetivo
Bienvenid@s a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Instrucciones:
- La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacerla de forma individual.
- La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
- El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
- Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
- No serán revisadas tareas desarrolladas en Python.
- Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
- Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.
Referencias:
Slides de las clases:
Videos de las clases:
- Introduction to Bayesian Inference: video1 video2 video3 video4 video5 video6 video7
- Summarizing the Posterior: video1 video2 video3 video4
Documentación:
Primera Parte: Preguntas Teóricas
A continuación, se presentaran diferentes preguntas que abordan las temáticas vistas en clases. Por favor responda cada una de estas preguntas de forma breve, no más de 4 o 5 lineas.
Pregunta 1:
Explique cual es la diferencia fundamental entre la estadística bayesiana y frecuentista.
El enfoque frecuentista requiere que todas las probabilidades sean definidas con respecto a frecuencias límite, es decir, sobre grandes muestras. Por tanto, en el enfoque frecuentista se trata como variable aleatoria sólo a cosas que se puedan medir. Los parámetros y modelos no pueden tener distribuciones de probabilidad, solamente los objetos medibles pueden tener una distribución de probabilidad. La distribución de estas medidas corresponde a la sampling distribution y sobre ella se puede realizar inferencia desde el punto de vista de la incertidumbre. Sin embargo, en la práctica casi siempre se trabaja con sólo una muestra. La escuela bayesiana plantear que, en realidad, una probabilidad es una creencia que puede estar basada en un supuesto como también en alguna percepción subjetiva. Por este motivo, en el enfoque bayesiano se relaja el supuesto de probabilidad por lo que es posible tratar a cosas no medibles como variables aleatorias. En el enfoque bayesiano se tiene una creencia inicial acerca de un parámetro, se obtienen datos en torno al parámetro, y en base a esa evidencia obtenida es posible actualizar la creencia inicial. El análisis de datos bayesianos no es más que contar los números en las formas en las que la data podría haber ocurrido de acuerdo a nuestros supuestos.
Pregunta 2:
Discuta la siguiente afirmación La inferencia bayesiana permite fácilmente utilizar distintos tipos de información.
La frase hace referencia a que a veces la información prior y los datos pueden ser de una naturaleza distinta. La información prior puede venir de un conocimiento que se codifica como un conteo compatible con los datos y puede ser actualizado de la misma forma, combinando el conocimiento prior subjetivo con los datos.
Pregunta 3:
Explique la diferencia entre prior probability y posterior probability
La prior probability corresponde al conocimiento a priori, a la creencia inicial que se tiene acerca de la distribución de los datos. Es posible mirar los datos y actualizar aquel conocimiento a priori. La posterior probability corresponde dicha actualización sobre la distribución de los datos, es decir, la posterior probability es conocimiento posterior que se tiene de la distribución de los datos dado el conocimento a prior (prior probability) y el conocimiento empiríco (los datos).
Pregunta 4:
El estadista bayesiano “Bruno Finetti” menciona la siguiente frase en su libro de probabilidad: La probabilidad no existe. Lo que en verdad quizo decir es que la probabilidad es un método para describir incertidumbre en un observador con conocimiento limitado. Discuta esta información utilizando el ejemplo del lanzamiento del globo terraqueo visto en clases. ¿Que significa decir “la probabilidad de que sea agua es 0.7”?
La frase se refiere a que la probabilidad objetiva como tal no existe, puesto que, según el enfoque bayesiano, una probabilidad respecto a algún fenómeno no es un valor puntual sino es una distribución generada a través de información a priori (supuestos sobre la distribución de los datos) e información empírica (un conjunto específico de datos obtenidos a través de un experimento). Respecto al ejemplo del globo terráqueo visto en clases, se tiene que decir que la probabilidad de que sea agua es 0.7 se refiere a que, dado los experimentos realizados, es decir, dado los datos, el valor mas probable para la proporción de agua con respecto al total de superficie de la tierra es \(0.7\). Sin embargo, dicha proporción es desconocida y, por tanto, el valor \(0.7\) es solo el valor mas probable en una distribución de valores.
la probabilidad se trata de una sensación psicológica de un observador con conocimiento limitado, ya que opera directamente sobre el grado psicológico de confianza de un individuo respecto a una cierta suposición. Ejemplo globo terráqueo Decir que la probabilidad de que sea agua es 0.7 significa que ..-
Pregunta 5:
¿ Que ventaja entrega que la distribución de la posterior este en la misma familia de distribución de probabilidad que la del prior?. De un ejemplo de alguna distribución que posee este comportamiento.
La ventaja que la distribución de la posterior esté en la misma familia de distribución de probabilidad que la del prior es que, en dicho contexto, se puede calcular analíticamente la distribución de probabilidad de la posterior cuando la posterior es la conjugada de la likelihood elegida. Por ejemplo, la distribución Beta es la conjugada de la distribución Binomial y, por tanto, la distribución de la posterior (Beta) está en la misma familia de distribución de probabilidad que la del prior cuando se utiliza una likelihood (Binomial) y, en dicho escenario, se puede calcular analíticamente la distribución posterior (Beta).
Pregunta 6:
Señale y explique los pasos de la grid approximation para obtener la posterior y responda las siguientes preguntas:
- ¿Cual de los pasos señalados nos permite obtener una distribución de la posterior mas precisa?.
- Cuales son las limitaciones de la grid approximation.
Pasos de la Grid Aproximation:
En primer lugar, se define la grilla, es decir, se debe decidir cuál será en rango de valores posibles para el parámetro a estimar y cuantos puntos (posibles valores del parámetro) se usarán para construir la distribución de probabilidad (posterior) del parámetro.
En segundo lugar, se computan los valores del prior para cada valor del parámetro en la grilla. Luego se obtiene la likelihood en cada valor de los parámetros. Despues se procede a computar la posterior no estandarizada para cada valor del parámetro, multiplicando el prior por la verosimilitud. Finalmente se debe estandarizar la posterior, dividiendo cada valor por la suma de todos los valores.
La elección de una buena distribución prior permite generar una distribución de la posterior mas precisa.
Una de las limitaciones de la grid aproximation es que escala de mala manera a medida que se aumenta el número de parámetros.
Pregunta 7:
¿ Por qué es necesario aprender a trabajar con muestras de la posterior?.
Es necesario aprender a trabajar con muestras de la posterior ya que lo que se quiere hacer es un análisis general del parámetro a estudiar. Cuando se procede a realizar el cálculo de una posterior se ocupa solo la información asociada a un dataset. La información de un dataset es como la información de solo un experimento y, en ese sentido, trabajar con muestras de la prior permite realizar más experimentos para obtener información y un análisis más general sobre el parámetro estudiado.
Pregunta 8:
Señale si las siguientes afirmaciones son verdaderas o falsas, en el caso que sean falsas justifique su respuesta:
Los point estimates de la posterior no entregan información relevante en un estudio.
Un intervalo de confianza es un intervalo dentro del cual un valor de parámetro no observado cae con una probabilidad particular, mientras que un intervalo de credibilidad es un rango de valores en el que se estima que estará cierto valor desconocido.
La principal ventaja de HPDI frente a un intervalo de credibilidad es que si la posterior no distribuye de forma normal, el HPDI será capaz de detectar los puntos de interés, mientras que un intervalo de credibilidad lo ignoraría al asumir simetría.
Depende: Cuando se tiene, por ejemplo, una posterior que distribuye de forma normal, calcular el MAP entrega información relevante ya que el MAP es la media de dicha distribución normal. Sin embargo, por ejemplo, para posteriores con distribuciones asimétricas o uniformes, los points estimates no entregan información relevante sobre dichas distribución.
Verdadero
Verdadero: Aunque es un poco más general la ventaja de HDPI ya que si la posterior no distribuye de forma simétrica (una distribución normal es un un caso de una distribución simétrica pero existen otras distribuciones simétricas como la distribución binomial, distribución t-student, distribución uniforme, etcétera), entonces HDPI es mejor que los intervalos de credibilidad.
Pregunta 9
Suponga que tiene dos especies de pandas. Cada una de las especies es igual de común y es imposible distinguirlas físicamente. Una de las diferencias entre las especies es el tamaño de sus familias. Si denotamos por \(\theta\) a la especie del panda se tiene que, cuando la especie es \(\theta = 1\) tiene dos bebes un \(10\%\) de las veces mientras que la especie \(\theta = 2\) tiene dos bebes un \(20\%\) de las veces, mientras que el resto de veces ambas especies tienen un solo bebe.
Suponga que usted esta intentando determinar la especie de un panda que que tiene como registro de nacimientos al conjunto \(D\), considere quela especie de un panda que acaba de dar a luz a dos bebes, es decir \(D = \{\text{dos bebes}\}\)
¿Cual es la probabilidad de que pertenezca a la especie \(1\)?
Suponga ahora que el mismo panda acaba de dar a luz y esta vez es solo un bebe. Calcule la probabilidad posterior de que el panda sea de especie \(1\). ¿Que cambia con el calculo anterior?
Suponga que le ofrecen hacer un test genético a su panda, como suele ser común con los test no es perfecto y le mencionan las siguientes características:
- La probabilidad de que correctamente idenfitique a la especie \(1\) es de \(0.8\)
- La probabilidad de que correctamente identifique a la especie \(2\) es de \(0.65\)
Se administra el test y se obtiene un resultado positivo a la especie \(1\). Sin utilizar la información en \(D\) calcule la probabilidad posterior de que su panda sea de la especie \(1\). Repita sus cálculos utilizando la información recopilada en \(D\). ¿En que varían sus resultados?
Se quiere calcula la siguiente probabilidad: \[ P(\theta=1 | D) = \frac{P(D | \theta = 1) \cdot P(\theta=1)}{P(D)} \] donde \(D = \{\text{dos bebes}\}\)
Por enunciado, se tiene las siguientre probabilidades: \[P(D | \theta = 1) = 0.1 \] \[ P(\theta=1) = 0.5 \] \[ P(D) = \sum_{\theta} P(D | \theta) \cdot P(\theta) = 0.1\cdot 0.5 + 0.2 \cdot 0.5 \]
Reemplazando en la fórmula original se tiene: \[ P(\theta=1 | D) = \frac{0.1 \cdot 0.5} { 0.1 \cdot 0.5 + 0.2 \cdot 0.5} = \frac{1}{3} \]
Es decir, la probabilidad de que el panda pertenezca a la especie 1 dado que el panda tuvo 2 hijos es \(\frac{1}{3}\)
Ahora, se quiere calcula la misma probabilidad: \[P(\theta=1 | D) = \frac{P(D | \theta = 1) \cdot P(\theta=1)}{P(D)}\] pero con \(D = \{\text{dos bebes y luego un bebe}\}\)
Por enunciado y por multiplicación de eventos independiente, se tiene las siguientre probabilidades \[ P(D | \theta = 1) = 0.1 \cdot 0.9 \cdot 0.5 \] \[ P(\theta=1) = 0.5 \] \[ P(D) = \sum_{\theta} P(D | \theta) \cdot P(\theta) = 0.1 \cdot 0.9 \cdot 0.5 + 0.2 \cdot 0.8 \cdot 0.5\]
Reemplazando en la fórmula original se tiene que: \[ P(\theta=1 | D) = \frac{0.1 \cdot 0.9 \cdot 0.5} { 0.1 \cdot 0.9 \cdot 0.5 + 0.2 \cdot 0.8 \cdot 0.5} = 0.36 \]
Es decir, la probabilidad de que el panda pertenezca a la especie 1 dado que el panda tuvo 2 hijos y luego tuvo 1 hijo es \(0.36\)
Nótese que, al introducir un nuevo dato (el panda tuvo un hijo), la probabilidad de identificar correctamente la especie del panda cambia, ya que el cálculo va a variar debido a que se va a actualizar la información que se tiene sobre el panda.
Se quiere calcular la siguiente probabilidad: \[ P(\theta=1 | D) = \frac{P(D | \theta = 1) \cdot P(\theta=1)}{P(D)} \] donde \(D = \{\text{el test identifico al panda como especie 1 }\}\)
Por enunciado, se tiene las siguientre probabilidades: \[P(D | \theta = 1) = 0.8\] \[ P(\theta=1) = 0.5 \] \[ P(D) = \sum_{\theta} P(D | \theta) \cdot P(\theta) = 0.8\cdot 0.5 + (1-0.65) \cdot 0.5 \]
Reemplazando en la fórmula original se tiene: \[ P(\theta=1 | D) = \frac{0.8 \cdot 0.5} { 0.8 \cdot 0.5 + 0.35 \cdot 0.5} = 0.695 \]
Es decir, la probabilidad de que el panda pertenezca a la especie 1 dado que el test identifico al panda como de la especie 1 es \(0.695\)
Ahora, se quiere calcula la misma probabilidad: \[ P(\theta=1 | D) = \frac{P(D | \theta = 1) \cdot P(\theta=1)}{P(D)} \] pero con \(D = \{\text{dos bebes, luego un bebe y el test identifico al panda como especie}\}\)
Por enunciado y por multiplicación de eventos independiente, se tiene las siguientre probabilidades \[ P(D | \theta = 1) = 0.1 \cdot 0.9 \cdot 0.85 \cdot 0.5 \] \[ P(\theta=1) = 0.5 \] \[ P(D) = \sum_{\theta} P(D | \theta) \cdot P(\theta) = 0.1 \cdot 0.9 \cdot 0.85 \cdot 0.5 + 0.2 \cdot 0.8 \cdot 0.35 \cdot 0.5\]
Reemplazando en la fórmula original se tiene que: \[ P(\theta=1 | D) = \frac{0.1 \cdot 0.9 \cdot 0.85 \cdot 0.5} { 0.1 \cdot 0.9 \cdot 0.85 \cdot 0.5 + 0.2 \cdot 0.8 \cdot 0.35 \cdot 0.5} = 0.563\]
Es decir, la probabilidad de que el panda pertenezca a la especie 1 dado que el panda tuvo 2 hijos, luego tuvo 1 hijo y finalmente que el test identifico al panda como de la especie 1 es \(0.563\)
Nótese que, al agregar datos, la probabilidad disminuye.
Segunda Parte: Elaboración de Código
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:
# Manipulación de estructuras
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.1.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(dplyr)
library(tidyr)
library(purrr)
# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout# Manipulación de varios plots en una imagen.
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine# Análisis bayesiano
library(rethinking)## Loading required package: rstan## Loading required package: StanHeaders## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:tidyr':
##
## extract## Loading required package: parallel## rethinking (Version 2.13)##
## Attaching package: 'rethinking'## The following object is masked from 'package:purrr':
##
## map## The following object is masked from 'package:stats':
##
## rstudentSi no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:
install.packages("rethinking")En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:
install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')Preguntas: Introducción a Grid Approximation
Las primeras dos preguntas de esta tarea tienen como objetivo introducirlos en la inferencia Bayesiana utilizando la técnica Grid Approximation para obtener una aproximación de la posterior. Al finalizar los problemas ustedes deberán ser capaces de visualizar los efectos que tiene el prior en la posterior, saber cómo realizar una Grid Approximation y comprender como utilizar Percentile Interval (PI) en una posterior.
Pregunta 1:
Considere el dataset “moneda.csv” donde se encuentran los resultados de un experimento lanzando una moneda, el objetivo de esta pregunta es estudiar mediante técnicas de inferencia Bayesiana el valor de la probabilidad de que salga cara, representado por el valor \(1\). Puede usar la librerira rethinking durante toda esta pregunta (si lo desea).
dataMoneda <- read.csv("moneda.csv", header = TRUE)Programe el metodo grid approximation para distintos tamaños de experimento. ¿Como van cambiando las curvas posterior?
Repita el mismo análisis anterior pero utilizando el método de Laplace (no necesita programar el método, puede utilizar la libreria “rethinking”). ¿Como se comparan con los resultados anteriores?.
Grafique la densidad de la posterior y encuentre la proporción de los siguientes defined boundaries:
- \((0, 0.4)\)
- \((0.4,0.7)\)
- \((0.7,1)\)
¿Como puede interpretar los resultados?
- Calcule un intervalo de credibilidad al \(50\%\), \(75\%\) y \(95\%\). ¿Como se puede interpretar los resultados? ¿Cual podría ser un problema al usar intervalos de credibilidad?.
- Genere los intervalos HDPI para \(50\%\), \(75\%\) y \(95\%\), compárelos con los intervalos de credibilidad de la parte anterior. ¿En que se diferencian los intervalos de credibilidad con los HDPI?.
Respuesta
Inicialmente se programó el método grid_aprox que permite generar grid approximations para distintos tamaños de experimentos dado un dataset dado. En particular, la función genera un grid approximation para un tamaño de experimento en particular. Del mismo modo, se creó la función auxiliar fun_data que permite formatear la data para la función grid_aprox.
fun_data <- function(data, sample_size = NULL){
if(!is.null(sample_size)){
data = data[sample(length(data), sample_size, replace = TRUE)]
}
sum_data_1 = sum(data)
total = length(data)
sum_data_0 = total - sum_data_1
return(list(
sum_data_0=sum_data_0,
sum_data_1=sum_data_1,
total=total
))
}
grid_aprox <- function(data, grid_size=100){
sum_data_1 <- data$sum_data_1
total <- data$total
# define grid
p_grid <- seq(from=0 , to=1 , length.out=grid_size)
# define prior
prior <- rep(1 , grid_size)
# compute likelihood at each value in grid
likelihood <- dbinom(sum_data_1 , size=total, prob=p_grid )
# compute product of likelihood and prior
unstd.posterior <- likelihood * prior
# standardize the posterior, so it sums to 1
posterior <- unstd.posterior / sum(unstd.posterior)
return(list(x=p_grid, y=posterior))
}Luego, para graficar las curvas para distintos tamaño de experimentos se creó utilizó la función plots_aprox. Del mismo modo, se implementó la función auxiliar add_curve que facilita la implementación de la función plots_aprox.
add_curve <- function(plot, df_i, k, plot_p){
if(k=="10"){
plot <- plot + geom_line(data=df_i, aes(x=x, y=y, colour="10"))
if(plot_p){
plot <- plot + geom_point(data=df_i, aes(x=x, y=y, colour="10"))
}
}
else if(k=="50"){
plot <- plot + geom_line(data=df_i, aes(x=x, y=y, colour="50"))
if(plot_p){
plot <- plot + geom_point(data=df_i, aes(x=x, y=y, colour="50"))
}
}
else if(k=="100"){
plot <- plot + geom_line(data=df_i, aes(x=x, y=y, colour="100"))
if(plot_p){
plot <- plot + geom_point(data=df_i, aes(x=x, y=y, colour="100"))
}
}
else if(k=="300"){
plot <- plot + geom_line(data=df_i, aes(x=x, y=y, colour="300"))
if(plot_p){
plot <- plot + geom_point(data=df_i, aes(x=x, y=y, colour="300"))
}
}
else if(k=="500"){
plot <- plot + geom_line(data=df_i, aes(x=x, y=y, colour="500"))
if(plot_p){
plot <- plot + geom_point(data=df_i, aes(x=x, y=y, colour="500"))
}
}
else{
plot <- plot + geom_line(data=df_i, aes(x=x, y=y, colour="1000"))
if(plot_p){
plot <- plot + geom_point(data=df_i, aes(x=x, y=y, colour="1000"))
}
}
return(plot)
}
plots_aprox <- function(data, fun_aprox, plot_p=TRUE){
sample_sizes <- c(10, 50, 100, 300, 500, 1000)
plot <- ggplot()
for(i in sample_sizes){
format_data <- fun_data(data, i)
x_y <- fun_aprox(format_data)
df_i <- data.frame(x=x_y$x, y=x_y$y)
k = as.character(i)
plot <- add_curve(plot, df_i, k, plot_p)
}
plot <- plot + xlim(c(0, 1)) +
guides(fill=guide_legend(title="Sample size")) +
scale_colour_manual("Sample size",
breaks = sapply(sample_sizes, as.character),
values = c("red", "blue", "green", "yellow", "gray", "black")) +
labs(x="Probabilidad de obtener cara",
y="Probabilidad de obtener x probabilidad",
title =paste("Aproximación con método ", as.character(substitute(fun_aprox))))
return(plot)
}Ahora, dada la función grid_aprox, se graficaron las curvas posteriores para distintos tamaños de experimentos.
plots_aprox(dataMoneda$Resultado, grid_aprox)
Nótese que, mientras mayor la cantidad de experimentos, la distribución tiende al valor de probabilidad asociada a los datos, es decir, se tiene que las distribuciones con mayores números de experimentos se comportan similar a una distribución normal con \(\mu\) la probabilidad de obtener cara dado los datos y valores pequeños para \(\sigma\).
Con la función laplace_aprox se repitió el procedimento que se realizó para el caso de la grid approximation.
laplace_aprox <- function(data, sample_size){
sum_data_0 = data$sum_data_0
sum_data_1 = data$sum_data_1
results <- quap(
alist(
W ~ dbinom(W+L, p),
p ~ dunif(0,1)
),
data=list(W=sum_data_1, L=sum_data_0))
results <- precis(results)
pop_mean <- results$mean
pop_sd <- results$sd
x <- seq(-6, 6, length = 1000) * pop_sd + pop_mean
y <- dnorm(x, pop_mean, pop_sd)
y <- y /sum(y)
return(list(x=x, y=y))
}
plots_aprox(dataMoneda$Resultado, laplace_aprox, plot_p=FALSE)## Warning: Removed 474 row(s) containing missing values (geom_path).## Warning: Removed 58 row(s) containing missing values (geom_path).
Del mismo modo que para el caso de la grid approximation, se tiene que, mientras mayor la cantidad de experimentos, la distribución tiende al valor de probabilidad asociada a los datos, es decir, se tiene que las distribuciones con mayores números de experimentos se comportan similar a una distribución normal con \(\mu\) la probabilidad de obtener cara dado los datos y valores pequeños para \(\sigma\).
Para graficar la densidad de la posterior se tiene el siguiente bloque de cídigo.
format_data <- fun_data(dataMoneda$Resultado)
aprox <- grid_aprox(format_data, grid_size = 1000)
samples <- sample(aprox$x, prob=aprox$y, size=1e4, replace=TRUE)
ggplot(data=data.frame(x=samples)) +
geom_density(aes(x=x)) +
# xlim(c(0.2,0.8)) +
labs(x="Probabilidad de obtener cara",
y="Probabilidad de obtener x probabilidad",
title = paste("Real posterior"))
Para calcular las proporción de los defined boundaries dados el siguiente bloque de código permite realizar dicha acción.
bounds <- function(data, low, upp, grid_size=100){
format_data <- fun_data(data)
aprox <- grid_aprox(format_data, grid_size = grid_size)
samples <- sample(aprox$x, prob=aprox$y, size=1e4, replace=TRUE)
prop <- sum(samples>=low & samples<=upp)/1e4
return(prop)
}
bounds(dataMoneda$Resultado, 0, 0.4, 10)## [1] 0bounds(dataMoneda$Resultado, 0.4, 0.7)## [1] 1bounds(dataMoneda$Resultado, 0.7, 1)## [1] 0Nótese que tanto para el intervalo \([0, 0.4]\) como para el intervalo \([0.7, 1]\) la proporción calculada es \(0\). Lo anterior tiene sentido con el gráfico de densidad de la posterior ya que, como se observa en el gráfico, los datos están netamente concentrados en una pequeña vecindad cercana al valor \(0.56\), es decir, el área total bajo la curva se encuentra, aproximadamente, entre los valores \(0.51\) y \(0.60\).
Ahora, para el cálculo de los intervalos de credibilidad se tiene lo siguiente:
credible_intervals <- function(data, p, grid_size=100){
format_data <- fun_data(data)
aprox <- grid_aprox(format_data, grid_size = grid_size)
samples <- sample(aprox$x, prob=aprox$y, size=1e4, replace=TRUE)
interval <- PI(samples, prob=p)
return(interval)
}
credible_intervals(dataMoneda$Resultado, 0.5)## 25% 75%
## 0.5454545 0.5656566credible_intervals(dataMoneda$Resultado, 0.75)## 12% 88%
## 0.5353535 0.5656566credible_intervals(dataMoneda$Resultado, 0.95)## 3% 98%
## 0.5252525 0.5858586Por el otro lado, para el cálculo de los intervalos HDPI, se tiene el siguiente bloque de código
hdpi <- function(data, p, grid_size=100){
format_data <- fun_data(data)
aprox <- grid_aprox(format_data, grid_size = grid_size)
samples <- sample(aprox$x, prob=aprox$y, size=1e4, replace=TRUE)
interval <- HPDI(samples, prob=p)
return(interval)
}
hdpi(dataMoneda$Resultado, 0.5)## |0.5 0.5|
## 0.5353535 0.5555556hdpi(dataMoneda$Resultado, 0.75)## |0.75 0.75|
## 0.5353535 0.5656566hdpi(dataMoneda$Resultado, 0.95)## |0.95 0.95|
## 0.5252525 0.5757576Respecto a los intervalos de credibilidad, un problema de ocupar dichos intervalos es que, para distribuciones poco simétricas, no son óptimos porque son simétrico respecto a su centro, es decir, el peso de cada lado es equivalente. En ese sentido, los intervalos HDPI sirven para tanto para distribuciones simétricas como no simétricas.
En el caso para el dataset de experimentos de lanzar una moneda, la diferencia entre los intervaos de credibilidad y los intervalos HDPI es pequeña, es decir, para consultas del mismo intervalo se obtiene resultados similares. Lo anterior tiene sentido ya que, al observar el gráfico de densidad asociado a dicho dataset, se evidencia una curva bastante simétrica y, por tanto, ambos intervalos funcionan correctamente.
De todas formas, como la curva de densidad no es exactamente simétrica, existen pequeñas diferencias en los dos tipos de intervalos utilizados (pero diferencia poco significativas).
Pregunta 2: Grid Approximation y Posterior Predictive Distribution
El objetivo de esta pregunta es comprender el concepto de sample prediction visto en clases y realizar predicciones en base a una posterior.
Un conjunto de carteros aburridos de las mordidas de perros ha decidido realizar un catastro de mordidas recibidas por los empleados de su empresa en un periodo de dos meses, planeando en base a estos datos realizar inferencia bayesiana. Los datos de las mordidas estas datos por el dataset no+mordidas.csv, en donde cada fila representa las mordidas recibidas por diferentes carteros y las columnas señalan si fue mordido o no el cartero en los meses de estudio (notar que si fue mordido sera señalado con un 1, de lo contrario es señalado con un 0). Cabe señalar que un cartero no puede ser mordido mas de una vez al mes, ya que el damnificado recibe licencia por todos los días restantes del mes tras la mordida, reincorporándose el siguiente mes al trabajo.
df = read.csv("no+mordidas.csv")
head(df)En base a los datos, realice los siguientes puntos:
Realice una grid approximation para estimar la probabilidad que un cartero sea mordido, para esto junte los datos del mes 1 y 2 de estudio. Señale el máximo valor de la posterior.
Utilizando la posterior obtenida en el paso anterior, utilice rbinom para simular 10.000 réplicas de 500 registros de mordidas. Con esto, deberá obtener 10.000 números, cada uno de los cuales es un recuento de las mordidas obtenidas en el registro de datos. Compare la distribución del número de los carteros mordidos predichos con el número real de los datos (248 carteros mordidos de un total de 500 datos). ¿El modelo se ajusta bien a los datos? Es decir, ¿la distribución de las predicciones incluye la observación real como resultado central y probable?
Como se comento en el comienzo
bites_month1contiene las mordidas señaladas por un conjunto de personas en el primer mes. Haciendo uso debites_month2, obtenga la posterior de que una persona que fue mordida en el primer mes, sea mordida nuevamente en el segundo mes. Para esto, se recomienda comenzar buscando los carteros que fueron mordidos el primer mes y en base a estos generar una búsqueda indexada para obtener el número solicitado. Hecho esto, simule ese número carteros mordidos 10.000 veces. De los resultados obtenidos, compare el recuento de carteros mordidos con el recuento real. ¿Cómo se ve el modelo desde este punto de vista?
Repuesta
Inicialmente, se juntaron los datos del primer mes con el segundo mes.
m1 <- df$bites_month_1
m2 <- df$bites_month_2
total <- append(m1, m2)Luego, se generó una grid approximation y se grafico la curva asociada a dicha aproximación.
plot_grid <- function(data, n, title){
plot <- add_curve(ggplot(), data, n, TRUE)
plot <- plot + xlim(c(0, 1)) +
guides(fill=guide_legend(title="Sample size")) +
scale_colour_manual("Sample size",
breaks = c(as.character(n)),
values = c("black")) +
labs(x="Probabilidad de ser mordido",
y="Probabilidad de obtener x probabilidad",
title =paste(title, "\nGrid Approximation"))
return(plot)
}
format_data <- fun_data(total)
aprox <- grid_aprox(format_data)
df_mor <- data.frame(x=aprox$x, y=aprox$y)
plot_grid(df_mor, format_data$total, "Distribución de probabilidad de que un cartero sea mordido")
Del gráfico se puede observa que se tiene que el valor con la mayor probabilidad posterior (también conocido como máximo a posterior o MAP) es \(0,4950\), es decir, la probabilidad más probable de que un cartero sea mordido por un perro es del \(49,50 \%\).
Lo anteriormente mencionado se comprueba numéricamente con el siguiente bloque de código.
x <- df_mor$x
y <- df_mor$y
max_index <- which(y %in% max(y))
max_p <- x[max_index]
max_p ## [1] 0.4949495Ahora, utilizando la posterior obtenida en el paso anterior, se utilizará la función rbinom para simular 10.000 réplicas de 500 registros de mordidas. Para graficar dicha información se, implementó la función plot_his
plot_his <- function(y_binom, title){
y_binom <- y_binom
plot <- ggplot(data.frame(x=y_binom), aes(x=x)) +
geom_histogram(binwidth = 2) +
labs(x="Número de carteros mordidos",
y="Frecuencia",
title=title)
return(plot)
}Utilizando solo en valor MAP se obtuvo el siguiente histograma.
y_binom_map <- rbinom(1e5, size=500, prob=max_p)
plot_his(y_binom_map, "Histograma nuevas predicciones utilizando MAP")
Por el otro lado, utilizando la distribución posterior se obtuvo el siguiente histograma.
samples <- sample(x, prob=y, size=1e4, replace=TRUE)
y_binom_samples <- rbinom(1e5, size=500, prob=samples)
title ="Histograma nuevas predicciones utilizando la distribución aproximada
de la variable p (probabilidad de que un cartero sea mordido)"
plot_his(y_binom_samples, title)
Nótese que el modelo se ajusta bien a los datos ya que la distribución de las predicciones incluye la observación real como resultado central y probable. Se puede observar que el valor MAP calculado sobre la grid approximation es central en los histogramas calculados y es de los valores mas probables. Del mismo modo, la forma (distribución) del histograma es similar a la distribución asociada a la curva de densidad. Sin embargo, la curva de densidad está más centrada sobre el MAP, es decir, en la curva de densidad solo valores muy cercanos al MAP tienen probabilidad significativa. En cambio, el histograma le da un rango más amplio y probable al valor de la probabilidad. Lo anteriormente mencionado es bueno ya que permite no centrar los resultados y el análisis en un dataset en particular sino que permite, dada el dataset, generar resultados más generales y tener conclusiones más abiertas sobre la probabilidad (lo cual desde un enfoque bayesiano es correcto ya que la probabilidad no en un valor en sí mismo)
Ahora, para el caso de los cartores mordidos los dos meses, el siguiente bloque de código calcula la cantidad de dichos carteros.
m1_1 <- which(m1==1)
m2_1 <- m2[m1_1]
m1_m2_1 <- m2[m2_1]
m1_m2_1 <- append(m1_m2_1, rep(0, 250-length(m1_m2_1)))Con los datos obtenidos, se generó una grid approximation y se grafico la curva asociada a dicha aproximación.
format_data <- fun_data(m1_m2_1)
aprox <- grid_aprox(format_data)
df_mor_1_2 <- data.frame(x=aprox$x, y=aprox$y)
plot_grid(df_mor_1_2, format_data$total, "Distribución de probabilidad de que un cartero sea mordido\nlos dos meses")
Del gráfico se puede observa que se tiene que el valor con la mayor probabilidad posterior (también conocido como máximo a posterior o MAP) es \(0,2222\), es decir, la probabilidad más probable de que un cartero sea mordido por un perro los dos meses es del \(22,22\%\).
Lo anteriormente mencionado se comprueba numéricamente con el siguiente bloque de código.
x <- df_mor_1_2$x
y <- df_mor_1_2$y
max_index <- which(y %in% max(y))
max_p <- x[max_index]
max_p ## [1] 0.2222222Utilizando solo en valor MAP se obtuvo el siguiente histograma.
y_binom_map <- rbinom(1e5, size=500, prob=max_p)
plot_his(y_binom_map, "Histograma nuevas predicciones utilizando MAP")
Por el otro lado, utilizando la distribución posterior se obtuvo el siguiente histograma.
samples <- sample(x, prob=y, size=1e4, replace=TRUE)
y_binom_samples <- rbinom(1e5, size=500, prob=samples)
title ="Histograma nuevas predicciones utilizando la distribución aproximada
de la variable p (probabilidad de que un cartero sea mordido)"
plot_his(y_binom_samples, title)
Respecto a la curva de densidad y los histogramas, las conclusiones son similares para cuando se juntó los datos de mordidas del primer y segundo mes. En primer lugar, el modelo se ajusta bien a los datos ya que la distribución de las predicciones incluye la observación real como resultado central y probable. Se puede observar que el valor MAP calculado sobre la grid approximation es central en los histogramas calculados y es de los valores mas probables. Del mismo modo, la forma (distribución) del histograma es similar a la distribución asociada a la curva de densidad. Sin embargo, la curva de densidad está más centrada sobre el MAP, es decir, en la curva de densidad solo valores muy cercanos al MAP tienen probabilidad significativa. En cambio, el histograma le da un rango más amplio y probable al valor de la probabilidad. Lo anteriormente mencionado es bueno ya que permite no centrar los resultados y el análisis en un dataset en particular sino que permite, dada el dataset, generar resultados más generales y tener conclusiones más abiertas sobre la probabilidad (lo cual desde un enfoque bayesiano es correcto ya que la probabilidad no en un valor en sí mismo)
Pregunta 3: Inferencia Sobre Dos Parámetros
En esta pregunta se trabajara con el dataset “notas.csv” el cual contiene las notas históricas de un curso desconocido. Suponga que los datos vienen de una distribución \(\mathcal{N}(\mu,\sigma^2)\), el objetivo de la pregunta es estudiar el comportamiento de los datos y los posibles valores de \(\mu,\sigma\) mediante técnicas de inferencia Bayesiana.
Usted sabe un dato extra sobre la información, los valores de \(\sigma\) en la grilla se mueven en el intervalo \([0.5,1.5]\) y, además, tiene una fuerte creencia a que es mas probable encontrar la desviación estándar real entre \([0.5,1]\) que en \((1,1.5]\). De hecho, estudios señalan que la probabilidad de encontrar sigma en los valores \([0.5,1]\) es de 2/3, mientras que 1/3 para el resto de intervalos.
- Modifique el siguiente código que permite realizar una grid approximation para \(2\) parámetros. Proponga priors para \(\mu\) y \(\sigma\), justifique su elección.
# Leer información
data_notas <- read.csv("notas.csv")
grid_size <- 100
# Función para crear likelihood dado mu y sigma
grid_function <- function(mu,sigma){
return(mean(dnorm(data_notas$Notas, mean = mu, sd = sigma)))
# return(dnorm(mean(data_notas$Notas), mean = mu, sd = sigma))
}
# Valores de la grilla
grid_mu <- seq(from=1 , to=7 , length.out=grid_size)
grid_sigma <- seq(from=0.5, to=1.5, length.out=grid_size)
# Se crea la grilla 2d
data_grid <- expand_grid(grid_mu,grid_sigma)
# Se guarda la likelihod
data_grid$likelihood <- map2(data_grid$grid_mu,data_grid$grid_sigma, grid_function)
# Se transforma el forma de map2 a una columna
data_grid <- unnest(data_grid,cols = c("likelihood"))
# Valores de los priors
prior_mu <- rep(1 , grid_size)
prior_sigma <- append(rep(2/3, grid_size/2), rep(1/3, grid_size/2))
# Se crea la grilla 2d de priors
prior <- expand_grid(prior_mu, prior_sigma)
# Se calculan los valores del prior
data_grid$prior <- map2(prior$prior_mu, prior$prior_sigma, prod)
data_grid <- unnest(data_grid,cols = c("prior"))
# Se calcula el posterior
data_grid$unstd_posterior <- data_grid$likelihood * data_grid$prior
# Se estandariza el posterior
data_grid$posterior <- data_grid$unstd_posterior/sum(data_grid$unstd_posterior)
# Se ajustan los valores de la posterior para que no sean valores tan pequñeos
data_grid$posterior <- (data_grid$posterior - min(data_grid$posterior))/(max(data_grid$posterior)-min(data_grid$posterior))- Tras haber ejecutado el código de la parte anterior ejecute el siguiente, ¿Que puede decir de los valores de la distribución? Se recomienda hacer modificaciones en el código para realizar un mejor análisis y estudio.
# Punto de referencia
# Se recomienda cambiar estos valores por unos adecuados que le permitan estudiar
# Los valores de la distribución de mejor manera
valor_x <- mean(data_notas$Notas)
valor_y <- sd(data_notas$Notas)
# Grafico
punto_comparacion <- tibble(x = valor_x, y = valor_y)
plt <- data_grid %>%
ggplot(aes(x = grid_mu, y = grid_sigma)) +
geom_raster(aes(fill = posterior),
interpolate = T
)+
geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
geom_label(
data = punto_comparacion, aes(x, y),
label = "Punto Comparación",
fill = "green",
color = "black",
nudge_y = 0, # Este parametro desplaza la caja por el eje y
nudge_x = 1 # Este parametro desplaza la caja por el eje x
)+
scale_fill_viridis_c() +
labs(
title = "Posterior para Mean y Standard Deviation",
x = expression(mu ["Mean"]),
y = expression(sigma ["Standar Deviation"])
) +
theme(panel.grid = element_blank())
plt
- A continuación se presenta un código que permite realizara la distribución dada por sampling from grid approximation ¿Para que sirve este proceso? ¿Que puede deducir del gráfico?
# Codificamos los datos
x <- 1:length(data_grid$posterior)
# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)
# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]
# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)
# Realizamos las densidades
dens(df)
1 - Al visualizar el gráfico, se puede decir que los valores mas posibles (color amarillo) de \(\mu\) se encuentran cercanos al valor \(6\) lo cual tiene sentido con respecto al promedio de los datos de las notas. Del mismo modo, para \(\sigma\) se puede realizar una afirmación similar: los valores mas posibles (color amarillo) de \(\sigma\) se encuentran cercanos al valor \(0.5\) lo cual tiene sentido con respecto a a desviación estándar de los datos de las notas. Asimismo, el efecto de las prior establecidas (tanto para \(\mu\) como para \(\sigma\)) se evidencia en el gráfico.
2- Este proceso de obtener muestras a partir de la grid aproximation permite obtener la posterior, con el objetivo de poder encontrar las distribuciones de la media y de la desviación estándar de las notas.
A work by CC6104